Explicit symplectic integrators for solving non-separable Hamiltonians. 
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I. INTRODUCTION 
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By exploiting the error functions of explicit symplectic integrators for solving separable Hamilto- 
nians, I show that it is possible to develop explicit, time-reversible symplectic integrators for solving 
, non-separable Hamiltonians of the product form. The algorithms are unusual in that they are of 

' fractional orders. 

o 
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^0 ' Symplectic integrators^ ^ are the methods of choice for solving diverse physical problems in classical^'^ ^, 
quantum*^^^, and statistical^^"^^ mechanics. For separable Hamiltonians, the problem is well understood and many 
p^j ■ explicit integrators are available^"''. However, for non-separable Hamiltonians, only implicit algorithms are known^"**. 
It is generally believed that no explicit algorithms can be developed for solving non-separable Hamiltonians'^'^. In this 
work, I show that this is not the case. Explicit, time-reversible algorithms can be developed to solve a selected class 
of non-separable Hamiltonians. The idea is to model non-separable Hamiltonians by the error terms of explicit algo- 
I rithms when solving separable Hamiltonians. By a suitable choice of factorization (or split) coefficients, the explicit 
O-i algorithm can be made to solve error Hamiltonians which are generally non-separable. 

^ , In the usual study of symplectic integratiors, one seeks to eliminate error terms in order to produce higher order 
Q ' algorithms. These error terms are therefore not of direct interest and are rarely studied in their own right. In this 
O work, these error terms are the non-separable Hamiltonians we seek to solve. The method can solve non-separable 
c/3 Hamiltonians of the product form, (sum over repeated indices) 
O 

= r,(p)F.,(q)T,(p), (1) 

■ provided that 



T,(p) = At(p). (2) 



and 



(N 
> 

■ For one degree of freedom, given T'{p) and V"{q), T{p) and V{q) can always be obtained by integration. 

I In the next section we will briefly summarize essential aspects of symplectic integrators and their error functions, 

. followed by our explicit integrator for solving the above non-separable Hamiltonian. Higher order algorithms are 

Q ' discussed in Section IV. 

>>: 

rH II. SYMPLECTIC INTEGRATORS 

Given a dynamical variable W{qi,pi) and a Hamiltonian function H(qi,pi), the former is evolved by the later via 
the Poisson bracket, and therefore by the corresponding Lie operator-^° H associated with the function H{qi,pi), 



= f^A_^Aw = w, (4) 

V dpi dqi dqi dpi / 



via exponentiation. 



For a separable Hamiltonian, 



W{t + e) ^e'"W{t). (5) 



i7(q,p) =r(p) + l/(q). 



(6) 
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the corresponding Hamiltonian operator is also separable, H = T + V, with T and V given by 

ay 9 

{;¥} = - — — . (8) 
Their corresponding evolution operators e^^ and e^^ then shift Qi and pj forward in time via 

, ^ ef 9T 

qi{e) =6" Qi = q, + s—, 
dV 

Pi{e) =e^^Pi=Pi-e—. (9) 

Conventional symplectic integrators correspond to approximating the short time evolution operator e^^ in the product 
form 

N 

e{f+V) _ TT ^Uef^VicV 



e 

2 = 1 



[]c*-^V'^^, (10) 



resulting in an ordered sequence of displacements (9) which defines the resulting algorithm. Here, we will consider 
only time-reversible, symmetric factorization schemes such that either ti = and Vi = Vn-i+i, U+i = tjv-i+i, or 
vn = and Vi = vn-i, U = tjv-i+i- 

The product of operators in (10) can be combined by use of the Baker-Campbell-Hausdorff (BCH) formula to give 



JV 



JJ^UeT^v.eV ^^eH^^ (11) 



j=l 

where the approximate Hamiltonian operator Ha has the general form 

Ha = e^f + e^V + e^e^^, [ffV] 

+e''ey^y[VfV]+0{e^) (12) 

where e^, e^^, ej.^^, etc., are functions of {U} and {vi} and where condensed commutator brackets, [TTF] = 
[T, [T, V]], [TVTV] = [T, [V, [T, V]]], etc., are used. From the way Lie operators are defined via (4), one can convert 
operators back to functions^'^ via [T, V] ^ {V,T} = —{T, V}, yielding 

HA=e^T + e^V + e^e^^^ {TTV} 

+e''e,^,{VTV} + 0{e^), (13) 

where again, condensed Poisson brackets, {TTV} = {T, {T, V}}, etc., are used. For a separable Hamiltonian of the 
form (6), we have 

dT dV 

By choosing {ti} and {vi} such that 

e^=ey = 0, (16) 
and either Cy^y = 0, or e^^y = 0, the algorithm would then be solving the non-separable Hamiltonian, either 

Httv^T,V,,T, or HvvT = V,T,,Vj. (17) 
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III. SOLVING NON-SEPARABLE HAMILTONIANS 

The following factorization scheme gives, 

T(e) = e'^''^^e'^*^^e^''^^e'^*i'^c^'""^c^*i*c^^i^c'^*^*e^''^^ 

= exp{s^[ffv]+e^E5+s'^E7 + e'^E9---), (18) 

with vo = — 2(fi + V2), t\ = —t2, V2 = —Vi/2 and Vi = l/t^- There is one free parameter t2 that one can choose to 
minimize the resulting error, but not be set to zero. As examplfied by (14) and (15), for a separable Hamiltonian 
H = T + V , higher order brackets of the form {T, Q}, {F, Q} have opposite signs. Thus one should choose algorithms 
with exQ = evQ to maximize error cancellations^^. This is the basis for symplectic corrector^^ or processed^^'^^ 
algorithms. The choice of t2 = —6^/^ « —1.82 forces bttttv = evTTTV and would be a good starting value. The 
RHS of (18) is the evolution operator for the non-separable Hamiltonian Httv with time step At = e'"* and leading 
error terms 0{e^). Thus the parameter e used by the integrator is £ = ^^At. Since = At^/^, the basic algorithm 
(18) in terms of At reads, 

T{At) = exp At {[ffv] + Ae'^E^ + At'^/^E^ + At^/^^g • • •) . (19) 

The order of the algorithm T{At) (the leading error in the Hamiltonian) is therefore only 2/3. We will discuss this 
and higher order algorithms in the next section. 

By interchange T F everywhere, but keeping the coefficents intact, the RHS of (18) goes over to 

^e^[ffV]^^e^[VVf]^ (20) 

and the basica algorithm T{At) solves the non-separable Hamiltonian HyvT- In both cases, the final force or velocity 
can be re-used at the start of the next iteration. Thus both algorithms require four-force and four- velocity evaluations. 
For one degree of freedom, any Hamiltonian of the form 

H = f{p)g{q) (21) 

can be solved. To test the algorithm, we solve the non-separable Hamiltonian 

HTTV = {l + ^)\l + q% (22) 

where the phase trajectory is harmonic near the origin, but highly distorted at larger values of (p, q). The algorithm's 
separable Hamiltonian is 

H=P+\p' + \q^ + ^q''. (23) 

In Fig.l we compare the phase trajectories produced by algorithm (18) with exact trajectories deduced from (22). 
We set t2 = —2 and use a relatively large value of At = 0.005 so that discrepances can be seen. The four trajectories 
are started at po = and g'o = 0.5, 1.0, 1.5, and 2.0 respectively. The error is largest at the positive maximum oip 
and next largest at the negative maximum of p. In each case, the error can be further reduced by making t2 more 
negative than -2. We did not bother with this refinement here, but this will be important in the 2D case discussed 
below. 

We will demonstrate that T{At) indeed converges as At^^^ in the next section. 
For more than one degree of freedom, the generalization of (21) to 

H = Y,fi{Pi)gi{qi) (24) 

i 

can always be solved. However, it is more interesting to generalize (23) to N-dimension by reinterpreting p and q as 
radial coordinates: p = ^/^~Pi, q = if - radial potential V{q), 

Vij = —6ij + (V" - — ) Mj , (25) 
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where here q is the unit vector. Thus the non-separable Hamiltonian Httv corresponding to the radial Hamiltonian 

(23) is 

Httv = (1 + y)' [l + + ^g'(p • q)'-] (26) 

This can again be solved by our explicit integrator (18). In two-dimension, most trajectories arc not closed and are 
likely to be chaotic. However, for some special initial configurations, a rich variety of closed orbits can be found. Fig. 2 
shows a sample of three such closed orbits. For this calculation, since the order of the algorithm is only 2/3, reducing 
the step size is not efficient in achieving higher accuaracy. Instead, we find that the error can be substantially reduced 
by changing t2 to w —3. For the circle, triangle and the twisted orbits of Fig. 3, the step sizes used were. At = 0.0012, 
0.001, and 0.0005 respectively. 

Finally, the standard kinetic energy term 

T{v) = \piPi (27) 



produces 



Httv = {TTV} = piVijpj , (28) 
HvTV = {VTV} = -ViVi, (29) 



and only Httv is non-separable. Here, Vij can be viewed as a position-dependent inverse mass matrix. This work 
shows that if Vij can be derived from a potential function V{q), then this non-scparablc Hamiltonian can also be 
solved by our explicit algorithm. Also, by itself, this quadratic Hamiltonian does not possess closed orbits for most 
V{q), thus explaining why this error term would disrupt closed orbit of the original Hamiltonian at large e. 

IV. HIGHER ORDER ALGORITHMS 

In the previous section, we have shown that the primitive algorithm T(Ai) does work and reproduces the correct 
phase trajectory. However, its 2/3-order convergence is very poor and requires extremely small At to produce accurate 
results. To demonstrate its fractional order convergence, we return to the one-dimensional case (22) and integrate 
from t = 0, po = 0, qo = 2 to t = = 0.385841, p{t) = —1.569196, q{t) = 0, corresponding to a quarter, clockwise 
rotation of the outermost phase trajectory of Fig.l. In Fig. 3, the relative error of the Hamiltonian (22) at t = T1/4 
is plotted as a function of At. The error of T{At) can be perfectly fitted with the power law — 2At^/'^, but due to 
this fractional power, the convergence at small At is very poor. Fortunately, the error structure (19) of T{At) allows 
simple ways of generating higher order symplectic algorithms. The triplet-construction of Creutz and Gocksch^^ and 
Yoshida^^ can produce arbitrary high order algorithms such as the following 4/3rd order algorithm 

--(-)-(^)-(-|^)-(^) (30. 
with s = 2^/^ and the following 6/3rd=2nd-order algorithm 

with s = 2^1'^ . As can be seen in Fig. 3, these higher order symplectic algorithms arc orders of magnitude better 
than the basic algorithm T{At). The disadvantage of the triplet-construction is that the computational efi^ort triples 
in going from order fc/3 to (fc + 2)/3. For example, the second-order algorithm 72(At) requires three evaluations of 
74/3 (Af), or nine evaluations of T(Af). Alternatively, arbitrary high order algorithms can also be obtained via the 
Multi-Product Expansion(MPE)^^, with only quadratically growing computational efforts. For example, by replacing 
kf — *■ fc^^^ in^^, one obtains 

^4r3""(A*) = ^n^t) + ) (32) 
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Tf'-^iAt) = - — ^ — TT^n^t) + 77, — 7^ — 

^ / 1 rj, On\/1n On.\ ^ ' (Or?, 1n\/On On\ 



(1" - 2")(1» - 3") ' ' (2" - 1»)(2» - 3") \2 



^ (in \n\!ni rm\^^ \ Q ) ('^"^) 



(3" - 1")(3" - 2") V 3 

with n = 2/3 in both cases. Here, T^ ^^^ {I\t) only requires six evaluations of T(A<). The disadvantage of MPE is 
that it is no longer symplcctic, but is like Rungc-Kutta-Nystrom type algorithms. However, as shown in Fig. 3, their 
energy error can be much smaller than the triplet symplectic algorithms. 



V. CONCLUDING SUMMARY 



In this work, we have shown that explicit symplectic integrators can be devised to solve a selected class of non- 
separable Hamiltonians. Any non-separable Hamiltonian which can be modelled by the error terms of an explicit 
integrator can be solved by the same integrator with changed splitting coefficients. The initial explicit algorithm is 
only of fractional order At^/^, but higher order algorithms can be easily obtained by use of the triplet construction 
or the multi-product expansion. 
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FIG. 1: The phase trajectories of the non-separable Hamiltonian (22). The computed phase points (stars) are compared with 
exact trajectories (hnes). The initial values are po = and qo = 0.5, 1.0, 1.5 and 2.0, corresponding to energy values of 1.25, 
2.0, 3.25 and 5.0 respectively. 
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FIG. 2: Some two dimensional orbits of the non-separable Hamiltonian (26). Most trajectory are not closed and only very 
special initial conditions can result in closed orbits. The initial conditions {qi,q2,Pi,P2) that produce the circle, the triangle 
and the twisted orbits are respectively, (0.8,0,0,0.425), (0.99,0,0,0.789) and (2.5,0,0,0.1884). 
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FIG. 3: The fractional power convergence of various explicit algorithms. The relative energy error is evaluated at the first 
quarter period t = 0.385841, for the outermost trajectory of Fig.l. The solid circles denote results of symplectic algorithms 
(18), (30) and (31). The hallow circles give results of MPE algorithms (32) and (33). The lines are fitted cuves of the form 
cAr, with n = 2/3,4/3, or 2 as indicated. 



